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It has recently been proposed that proteins embedded in lipidic bio-membranes can spontaneously 
self-organize into stable small clusters, or membrane nano-domains, due to the competition between 
short-range attractive and longer-range repulsive forces between proteins, specific to these systems. 
In this paper, we carry on our investigation, by Monte Carlo simulations, of different aspects of 
cluster phases of proteins in bio-membranes. First, we compare different long-range potentials 
(including notably three-body terms) to demonstrate that the existence of cluster phases should be 
quite generic. Furthermore, a real membrane contains hundreds of different protein species that 
are far from being randomly distributed in these nano-domains. We take this protein diversity 
into account by modulating protein-protein interaction potentials both at short and longer range. 
We confirm theoretical predictions in terms of biological cluster specialization by deciphering how 
clusters recruit only a few protein species. In this respect, we highlight that cluster phases can turn 
out to be an advantage at the biological level, for example by enhancing the cell response to external 
stimuli. 



1. INTRODUCTION 

Live cells contain a large amount of membranes play- 
ing a great variety of biological roles. Membranes are an 
essential component of Life. Their first function is to sep- 
arate the diverse cell organelles (e.g. the Golgi apparatus, 
mitochondria, the endoplasmic reticulum, the nucleus in 
animal cells) from each other. The plasma membrane, 
on which this paper primarily focusses, separates the in- 
terior of the cell from the environment. However, mem- 
branes are not passive barriers but they control in an 
active way the movement of molecules or information in 
and out the organelles or the cell, thus maintaining sub- 
stantial differences of composition between the cell and 
its environment. Plasma membranes also play an impor- 
tant role in cell adhesion and motility, thus ensuring a 
large amount of biological functions. The principal con- 
stituents of plasma membranes are lipids and proteins, 
accounting each for typically 50 % of the membrane mass, 
and constituting a two-dimensional complex fluid where 
interactions are weak, typically on the order of the ther- 
mal energy k^T [TJ [3J. Adopting a coarse-grained point 
of view, lipids can be considered as an underlying fluid 
(or solvent) in which included proteins evolve. Effec- 
tive protein-protein interactions mediated by this fluid 
emerge, because an inclusion creates a perturbation of 
the membrane that in turn influences neighbor inclusions. 

Membrane organization is a topical issue in cell biol- 
ogy. Understanding how hundreds of different protein 
and lipid species organize in a same membrane to per- 
form as many biological functions ^ [3] remains a chal- 
lenge to which physicists have recently paid attention, by 
applying concepts of soft matter to the biophysics of bio- 
membranes [5J H] . Far from the pioneer mosaic model of 
Singer and Nicolson [5], it is now understood that lipids 
and proteins are non-homogeneously distributed in mem- 
branes. They are subtly organized in small domains, or 



compartments, or clusters [6], of highly variable composi- 
tion, the length-scale of which ranges from a few nanome- 
ters to microns. Heterogeneity is even thought to be a 
key ingredient for biological functions by concentrating 
a few protein species in a same nano-domain. Indeed, 
it facilitates the encounter of different proteins species 
that must interact to perform together a given biological 
task [7 -10J. When in a same membrane domain, they 
encounter much more easily than if they were randomly 
diffusing on the membrane surface Concentrating 
identical receptors in a same domain also facilitates the 
response to small inputs by reducing noise |12j . Group- 
ing together structure proteins is also essential in the 
context of cell adhesion [T31 [2] . A wide range of exper- 
imental techniques demonstrate that co-localization of a 
few different protein species in the same membrane nano- 
domains is a common feature (see, e.g., [51 IT21 IT51 [PorBS] ) 
and as stated by T. Lang and S.O. Rozzoli, clusters are 
"more than a pretty picture" [H]. 

However, the precise way proteins organize and the 
physicochemical mechanisms responsible for this organi- 
zation are currently a matter of controversy [321 [23l l26l - 
150"] . In particular, a consensus has not been reached yet 
on the reason why domains remain sub-micrometric. The 
ideas developed in this paper propose elements of answers 
to the two following questions: 1) What causes compart- 
mentalization of membrane proteins in sub-micrometric 
domains in live cells? 2) How can compartmentaliza- 
tion mechanisms account for segregation of membrane 
proteins in specialized, heterogeneous sub-micrometric 
structures where they more easily encounter to perform 
biological functions? 

Motivated by this biophysical context, statistical me- 
chanics arguments have recently been proposed to ac- 
count for the existence of protein nano- or micro-domains 
in a membrane at equilibrium [H J l24 l [27 1 13TH33] . They 
rely on the following mechanism: attraction at short 
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range (~ 1 nm) between proteins favors condensation 
of membrane proteins in a dense phase, but some weaker 
repulsion at longer range (~ 10 nm) prevents a com- 
plete phase separation because once clusters have nucle- 
ated, sufficiently large clusters repel each other due to 
additivity of the repulsion between their proteins. Alter- 
natively, too large clusters are unstable because of the 
long-range repulsion. The resulting phase at equilibrium 
is called a "cluster phase" [32 , by analogy with simi- 
lar phases in soft condensed matter [34H36], However, 
we will show that real protein-protein interaction poten- 
tials are very complex and potentially very diverse. The 
present work addresses the question of the robustness of 
the cluster phases with respect to the potential shape 
and demonstrates that cluster phases should be generic 
in bio-membranes in spite of potential diversity. 

In addition, previous studies did not take into account 
protein diversity. Recent theoretical investigations have 
shown that modulating the short-range attraction in or- 
der to account for diversity could lead to composition 
heterogeneity of clusters because the ensuing energy gain 
is larger than the corresponding entropic cost [37 . The 
present work also addresses this question on numerical 
grounds and confirms these theoretical investigations (see 
also [EH], Supporting Material). 

The paper is thus organized as follows. We recall the 
biophysical context and we present our model in the two 
first sections. Then Section 4 is mainly devoted to the 
question of the robustness of cluster phases with respect 
to the potential shape, more precisely to its long-range 
repulsion term. In Section 5, we focus on the role of 
protein diversity through its effects on the short-range 
attraction. The last section is devoted to discussions and 
conclusions. 



2. BIOPHYSICAL CONTEXT 

Interactions between proteins embedded in bio- 
membranes are manifold and the literature on the topic 
is extensive. In addition to hard-core and electrostatic 
interactions, they feel several effective forces specific to 
these systems, because they are mediated by the mem- 
brane. As we shall see, the free energies involved in these 
interactions are on the order of the thermal energy k^T 
and thus play a role at equilibrium. 

If the membrane is seen as an elastic sheet of curvature 
elastic modulus K, it is locally deformed by the pres- 
ence of a protein, whether it be embedded or periph- 
eral, and the response to this constraint implies effec- 
tive long-range forces. First, thermal undulations of the 
elastic membrane in the dimension normal to the mem- 
brane plane are perturbed by the presence of proteins, 
in a way that depends on their separation r. The en- 
tropy, depending on the number of effectively accessible 
undulation modes, thus depends on separation, and an 
attractive Casimir-like force ensues, which is pairwise ad- 
ditive and decays as r -4 at large r [ESI HO]- In addition, 



when proteins are (up-down) asymmetric with respect to 
the membrane plane, which we assume to be the major- 
ity case for both peripheral and transmembrane proteins, 
they can be modeled as conical inclusions with a half- 
aperture angle (or contact angle) denoted by a. This 
angle models not only the crystallographic shape of the 
protein [13 02], but also its local interactions, in partic- 
ular of electrostatic nature, with the surrounding lipids. 
The two lipid monolayers of cell membranes are very dif- 
ferent in composition and in charge [1] , which also breaks 
the up-down symmetry, even for cylindrical inclusions. 
Such asymmetric inclusions feel an additional repulsion, 
also decaying as Kr~ A at leading order. This repulsion, 
proportional to the squares of the effective contact an- 
gles a, is due to the elastic deformation imposed to the 
elastic membrane [ESI SOI BEll45| . Note the important 
following point: even for differently oriented conical in- 
clusions with contact angles ot\ and «2 of opposite sign, 
this elastic contribution remains repulsive because it is 
proportionnal to a\ + a% at leading order. This repul- 
sion compensates the Casimir-like attraction as soon as 
a is typically larger than 5°. It is not pairwise additive 
anymore, since three-body terms, which can be attrac- 
tive, exist at the leading order r -4 [4"4ll46] . But we will 
demonstrate that they are not strong enough to counter- 
balance repulsive two-body forces at long range (see also 
the Supporting Information). As for n-body interactions 
terms with n > 3, they decay faster than r -4 [H] and 
will not be considered in the present work [45j . Further- 
more, proteins need not be isotropic inclusions, in this 
sense that their horizontal section can be better mod- 
eled by an ellipse (instead of a disk). In this case, the 
potential depends on the relative orientation of the el- 
lipses [miH], and decays as r~ 2 . But except for strong 
anisotropics, it is averaged over orientations because of 
the rapid protein rotational diffusion 41j, which signifi- 
cantly lowers the strength of the interaction and makes 
it decaying as r -4 as well. In case of strong anisotropy, 
the ensuing elastic interaction can even become attrac- 
tive |41j . In the present work, we only consider weakly 
or moderately anisotropic inclusions, such that the net 
elastic interaction is repulsive at long range. We dis- 
cuss strong anisotropy in the Supporting Information. 
Beyond anisotropy, protein shape can make interaction 
potentials even more complex [47] . 

Even though these interactions mediated by the mem- 
brane should decay algebraically at long range, they are 
likely to be partially screened for at least two reasons: 
1) Long wavelength membrane excitations are damped 
when the membrane is coupled to the much more rigid cy- 
toskeleton [JS]; 2) When weak, but non- vanishing mem- 
brane tension is taken into account, the repulsion is 
screened beyond distances of a few tens of nanometers 
and thus decays exponentially at long range [IS] (see 
also [31] [50] for related calculations). However, estimat- 
ing the values of membrane tensions in live cells, and 
therefore of the screening length, is a difficult experimen- 
tal task. It is thus relevant to consider both exponentially 
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and algebraically decaying potentials. 

At shorter range, proteins first experience a hard-core 
repulsion at contact. In this respect, protein radii a are 
variable, ranging from a fraction of nanometer for sin- 
gle membrane-spanning a-helices, to a couple of nanome- 
ters for large proteins (such as receptors, pumps or chan- 
nels) pp. As a first step, we consider mono-disperse dis- 
tributions of protein radii a in this work, with the typical 
radius a — 2 nm. The general case is briefly discussed in 
the Supporting Information. At contact or close to con- 
tact, interactions are of hydrophobic, or (screened) elec- 
trostatic, or hydrogen-bond origins, arising from apolar, 
or polar, or charged amino-acids near the protein surface. 
Because of the rotational averaging discussed above, we 
do not take into account the anisotropy of these interac- 
tions either. Van der Waals forces (in r -6 ) should also 
play a role at short distances. 

Two additional forces mediated by the lipidic mem- 
brane play a role when separations between protein sur- 
faces are typically less than one nanometer. The hy- 
drophobic core of transmembrane proteins is likely not to 
match the equilibrium membrane thickness, what is usu- 
ally called hydrophobic mismatch. The mismatch is said 
to be positive when the core is too thick, and negative 
in the converse case. Two proteins with mismatches of 
identical (resp. different) sign feel an attractive (resp. re- 
pulsive) short-range force due to the induced elastic vari- 
ations of the bilayer thickness p~5j l5lTf54] . Mismatch has 
been experimentally demonstrated to promote the ag- 
gregation of transmembrane proteins |55j . There also ex- 
ists many-body effects associated with hydrophobic mis- 
match [56] . They are not included in the present study 
because no simple effective expression is known for these 
forces. 

In addition, attractive depletion forces, due to the two- 
dimensional osmotic pressure of lipids on proteins, tend 
to bring them closer when they are about a nanometer 
apart 57J. Hydrophobic and depletion forces have been 
studied numerically by Molecular Dynamics [3SJ E3J 1551 - 
160] . Their range is nanometric and the involved energies 
are equal to a few k^,T. Attractive or repulsive hydropho- 
bic mismatch modulates this strength by a couple of k-gT 
and can consequently change the degree of aggregation 
(or oligomerization) of proteins (e.g. rhodopsins in |61|). 

Finally, cell membranes are constituted of a large vari- 
ety of different lipid species |62| and proteins are known 
to recruit in their neighborhood lipids for which they 
have a better affinity [3J. For two identical proteins, the 
closer they are, the more energetically favorable the con- 
figuration because there is an interface energy associated 
with this "wetting" phenomenon. Typical energies are 
also of order k^T [5TJH2!, 63 . These references propose a 
protein-driven mechanism for domain formation invoking 
such forces, but the limited domain size due to additional 
repulsive forces has not been discussed in this context. 

To conclude, protein interaction potentials have di- 
verse contributions, and writing a generic potential shape 
is a tedious task appealing to intensive Molecular Dynam- 



ics simulations, out of the scope of the present work. For 
this reason, it is useful to consider different realistic, typ- 
ical shapes, and to analyze the similarities and the differ- 
ences between the ensuing cluster phases, as already dis- 
cussed in reference [35] • We focus on two types of interac- 
tion potentials: potentials decaying algebraically as r~ 4 
at long range, corresponding to unscreened or partially 
screened repulsion; and potentials decaying exponentially 
at long range, corresponding to complete screening. We 
thus span a large variety of experimental contexts. We 
demonstrate that conclusions are qualitatively similar in 
both cases, thus proving that cluster phases of membrane 
proteins should be generic. 



MODEL 



Here and in the following, all 
units of k^,T, because it is the 
in this context. 

The pairwise potentials Viir 
comprise a hard-core repulsion 
i.e. Ui(r) — oo if r < 2a, an 
range and a weak repulsive one 
attractive part is chosen to be 
because the short-range forces 
rapidly by nature. 



energies are implicitly in 
relevant scale of energies 

) considered in this work 
at very short distances, 
attractive part at short 
at larger distances. The 
exponential in all cases, 
considered above decay 




8 10 12 14 16 18 20 

r (nm) 



FIG. 1: Examples of pair potentials U? used throughout this 
paper. Fro m le ft to right at small r: two algebraic potentials 
[equation fOJ] with K = 150 (black) and_50fc B T (red); and 

= 0.4 



two bi-exponential potentials [equation ( 3.1 1] with E m 
(blue) and O.lksT (green) and l/7r = 4 (blue) and 16 nm 
(green). E m [ n — —AkuT, and l/7a = 2 nm in all cases. In 
all cases, cluster co-exist with monomers at equilibrium, as 
illustrated by the four random snapshots (at density (f> — 0.1, 
same colors as the corresponding curves). 



As motivated above, we will study on the one hand 
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bi-exponential potentials of the form 



decaying as r 4 : 



^2 XP ( r ) = ~ £ a exp(-7 a r) + e r exp(-j r r), 



(3.1) 



where all parameters are positive. They contain an at- 
tractive and a repulsive term and "f r < 7 a (see exam- 
ples in figure [I]) . They correspond to complete screen- 
ing. Parameter values will be specialized in our nu- 
merical studies, in agreement with the biophysical con- 
text of the previous Section. As in reference |32j . we 
set l/7a = 2 nm, l/j r — 4 or 16 nm, and e a and e r 
are adapted to fit the desired energy at contact E m i n = 
C/2 Xp (2a) and the desired energy barrier, E max . For ex- 
ample, l/7 r = 4 nm, E nlln = —Ak^T and E max — 0.4fceT 
lead to e a — bb.lk^T and e r = 9.4fceT. The maximum 
value i?max is then reached at r = 9.8 nm. 

On the other hand, we will consider unscreened po- 
tentials, which still display an exponential short-range 
attraction but the repulsive part of which is slowly, alge- 
braically decaying: 



7- r pOW/ \ 



-e a exp(-7 Q r) + E S if r < r, 
if r > ri ' 



(3.2) 



where the energy shift E s enforces the continuity of U at 
rj. The existence of cluster phases for p — 2 has been 
explored in reference [35]. They also exist for p = 6 
(data not shown). Here we focus on p — 4 because 
the repulsion between isotropic inclusions in an elas- 
tic membrane displays such a long-range behavior, with 
C = S 2 [^f(a\ + a\) - £], where S = Tra 2 is the inclu- 
sion area, K is the membrane elastic modulus and a.\ t 2 
are contact angles [39j [40] . The first term in the squared 
brackets corresponds to pure elastic deformations of the 
membrane and the second one to the Casimir-like forces. 
The constant C is positive for sufficiently large values of 
K and aa.,2- The values of the contact angles a are diffi- 
cult to estimate experimentally. Here, we shall consider 
two typical values, a — or 10°, so that, with the val- 
ues of K considered below, C is always positive except 
if ot\ — «2 — 0. As above, we set 1/7 Q = 2 nm and we 
choose ri = 2a + l/-f a , l/-f a being the typical extend of 
the attraction beyond the hard core distance, 2a. Finally, 
e a is set so that E min has the desired value (see below). 
Examples are given in figure [T] 

As far as the elastic modulus K is concerned, the ex- 
perimental values found in the literature are very vari- 
able, because they depend strongly on the lipidic com- 
position, in particular the cholesterol concentration. The 
typical values of K range from 10 to 200/cb? 1 ([55HSg] 
and references therein). In this work, we consider the 
values K — 30, 50 and 150fceT, thus spanning the ex- 
perimental interval. Accordingly, if 01=0^2 = 10°, then 
£ max = 0.21, 0.40 and 1.34fc B T, respectively. 

In Section |4j we will also explore the role of long-range 
three-body forces mediated by the elastic membrane, also 



^pow = B 



E 

p<q<r 



cos(27 P ) cos(27 Q ) cos(27 r ) 



^pq^qr 



^pr^qr 



(3.3) 

where j p is the angle (r pq , v pr ), r pq = r q —r p , and so forth, 
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) [Hl[45]. When wished, U$ ow is 



and B 

added to the two-body contibutions of equation (13. 21). In 



this work, we only consider two- and three-body forces 
because they dominate at large separations. The more 
general case of many-body forces will be addressed in the 
Discussion. 

In order to take protein diversity into account, we 
also consider in section [5] the case were several protein 
species are present in the membrane. Forces then de- 
pend on the nature of the interacting proteins. On the 
one hand, we can play on the long-range terms by mak- 
ing e r in (3.1 1 or C in (3.2) depend on interacting protein 
species, through their contact angles a. On the other 
hand, we can modulate the short-range attraction by us- 
ing species-dependent values of e a and e r to modulate 
E min at fixed j a , 7r and E 

max ■ 

Note that so far, we have used the terminology "poten- 
tials or forces at long range" . However, our potentials are 
not long-ranged in the strict statistical mechanics sense 
of the word because, even though their range extends be- 
yond several tens of nanometers, the integral J rU{r)dr 
is convergent at large r. However, we shall keep the ter- 
minology "long-range" in the following, by opposition to 
"short-range" . 

We carry out Monte Carlo simulations in the canoni- 
cal ensemble as prescribed in reference [35] . More details 
can be found in the Supporting Information. Our systems 
contain TV > 100 particles that interact via the potentials 
discussed above, and evolve in a two-dimenional continu- 
ous medium of area A representing the lipidic "sea" (the 
solvent). Boundary conditions are periodic. The hard- 
core diameter is chosen as do = 2a — 4 nm, the typical 
diameter of a protein of average molecular weight [T]. 
Figures [T] [6] [9] and S2-3 display several simulation snap- 
shots for various parameter sets. The distributions of 
cluster sizes P{k) displayed below are defined as 



P(k) = {N k )/J2{N k ), 



(3.4) 



k=l 



where (Nk) is the measured average number of clusters of 
size k. The mean cluster size is defined as (k) = kP(k). 
It is also equal to the number of particles, N, divided by 
the number of clusters (including monomers) . 



4. LONG-RANGE POTENTIALS 

In this section, we focus on the role of long-range repul- 
sive potentials on cluster phases. As motivated above, we 
first compare two kinds of pairwise potentials: exponen- 



tially [equation (3.1)] and algebraically [equation (3.2)] 
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decaying ones. In a second time, we study the role of 



three-body forces at long range [equation (3.3)]. Finally, 
we focus on the implication of protein diversity as far as 
long-range forces are concerned. More precisely, we take 
into account the fact that a macroscopic fraction of up- 
down symmetric proteins might not be concerned by the 
elastic long-range repulsion. 



4.1. Comparison between two kinds of long-range 
pair potentials 

We first compare exponentially decaying forces and al- 



gebraically decaying ones (U(r) 



), associated with 



the repulsion due to the elastic deformations of the mem- 
brane. The N proteins are identical. Bi-exponential po- 
tentials have already been numerically studied in detail 
in references [321 EH [65] . The main feature of the ensu- 
ing cluster phases is that cluster-size distributions P(k) 
appear to be bimodal for broad ranges of parameters and 
concentrations: large, dense clusters co-exist with a gas of 
monomers. Figure [8] below will provide related distribu- 
tion examples. In addition, let us denote by 4> = Nd^/ A 
the (reduced) protein density. Above a limiting value </> c , 
a nice approximate proportionality regime, (k) ~ 4>/4> c : > 
can in general be observed on a few decades in this bi- 
exponential case. These observations have been given 
theroretical interpretations in references [331 ES] • 

For particles experiencing an algebraic repulsion at 
long range, the situation is somewhat different. Clus- 
ters still appear to be stable at equilibrium, since they 
exist independently of the initial configuration: an initial 
gas partially condenses into a cluster phase; converesly, 
a unique initial condensed droplet splits up into smaller 
clusters (see Supporting Information, figure S3). And 
when algebraic and bi-exponential potentials are globally 
comparable, as the red and blue potentials in figure [T] 
systems also look similar, as illustrated by the snapshots 
in the same figure (see also the insets of figure [2] where 
values of (k) are plotted). 

However, cluster-size distributions P(k) are not bi- 
modal any longer in the algebraic case, but they rather 
display a broad power-law-like regime with an expo- 
nent close to —2, until a maximum size, as illustrated 
in figure [2] (we have nevertheless observed that a bi- 
modal regime is restored for large bending moduli K and 
large (j> (see figure [7])). Thus can we speak of "cluster 
phases"? Defining unambiguously a cluster phase is not 
an easy task because it is not even clear that such phases 
are characterized by a true thermodynamic transition at 
4> c [371 ESI- And density fluctuations in gases can also 
lead to transient small multimers and to a value of (k) 
that can be slightly larger than 1. However, upon some 
approximations, it has been shown that density fluctu- 
ations in a gas phase, lead to a rapidly, exponentially 
decaying distribution |37j . In the present case, we rather 
observe a power-law-like regime instead, and maybe more 
importantly, we infer from the different data sets used in 



10" 



10" 



10" 



o 10- 



10 



10" 



g 10~ 1 



"D _2 
CD 10 



1 \ 

\ N 


10 
8 

j 6 


V . . O" ' 


\\ N 


' 4 
2 


C\ ( 


) 0.05 0.1 0.15 


0.2 


* \v 
V 


v v \ \ 
V \ \-> 






V \ x 

\\ \ 

\ \ 





10 

Cluster size 



100 



3.2 p 

2.8 - 

2,4 - 

2 - 



10" 



10 




3.6 3.7 3:8 3:9 4- -4.1 4.2 4.3 



10 

Cluster size 



100 



FIG. 2: Distributions of cluster sizes P(k) in the case of alge- 
braic long-range repulsion. The dot-dashed lines have slope 
-2, for comparison. Top: for three different sets of parame- 
ter values: K = 30fc B T, Emm = — 3.5fc B T and thus i5 m ax = 
0.21fc B T (full black lines); K = 50fc B T, Emm = -3.8fc B T 
and -Emax = 0.40fc B T (dashed red lines); and K = 150fc B T, 
Emm = -4fc B T and _E max = 1.34fc B T (dotted blue lines). And 
1/ja = 2 nm in all cases. From left to right at high cluster 
size, for each set of parameters, </> = 0.05 (N — 100), (j> = 0.1 
(N = 100) and = 0.2 (N = 200 except for K = 150k B T 
where N = 100). Log-log coordinates. Inset: Average cluster 
size (k) as a function of the density (j>. Same color and line 
styles as above. For comparison, we have also plotted in this 
inset (filled black circles) (k) for the bi-exponential potential 
(blue curve) of figure [l] Bottom: K = 50fc B T, N = 100 and 
4> = 0.05; From left to right at high cluster size: Emin = —3.7, 
-3.8, -3.9, -4, -4.1 and -4.2fc B T. Inset: Average cluster 
size (fc) as a function of -Emm- 



figure [2] that more than 25% (resp. 75%) of the proteins 
dwell in clusters containing more than k = 10 particles 
as soon as (k) > 3 (resp. (A;) > 7). Independently of 
a precise definition of a cluster phase as in [371 112] , this 
indicates that a macroscopic fraction of proteins live in 
assemblies, which is the biological mechanism we are pri- 
marily interested in. 

The proportionality regime is still being observed, even 
though on shorter concentration ranges (see figure [2j 



() 



top, inset, and figure [3j). At the highest density stud- 
ied, <f> = 0.2, average cluster sizes (k) range from 4 to 
8 for the parameter sets studied and clusters of size up 
to k — 20 are commonly observed. Note that lowering 
E m m at fixed K (i.e. increasing the attraction) should 
increase (fc) [55] . Unfortunately, due to the long-tail dis- 
tribution, increasing (k) requires increasing N, which is 
rapidly limiting in terms of computational cost. How- 
ever, we have studied the effect of varying |-E ro in| in the 
K = 50fc B T and = 0.05 case (figured bottom). The 
power-law-like regime extends as |-E m inJ grows and (k) 
increases concomitantly (see Inset). 
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FIG. 3: Mean cluster size (k) as a function of the density <j>. 
Comparison bewteen simulation data for algebraic long-range 
repulsion, with K — 150fcBl\ l/7a = 2 nm, E m i n — — Ak^T, 
and N = 100, and experimental measurements, inferred 
from [70] (see text). 



To conclude, this numerical study indicates that re- 
placing and exponentially decaying repulsion at long 
range by an algebraic one does not alter the existence of 
cluster phases for wide ranges of parameters, thus demon- 
strating that they should be robust with respect to the 
precise shape of interaction potential, provided that it 
displays both a short-range (~ 1 nm) attraction and a 
longer-range (~ 10 nm), weaker repulsion. 

To finish with, we compare our data to available exper- 
imental ones. Clusters of bacteriorhodopsin (BR) in pro- 
teoliposomes can be observed by freeze-fracture electron 
microscopy. This technique makes possible the determi- 
nation of cluster-sizes and of their distributions with a 
correct accuracy [7U] (see Table I and figure 4 in this ref- 
erence). We assume a BR radius of 1.7 nm [70] to convert 
numbers of BR per fim 2 to densities <f>. The algebraically 
decaying potential (3.2) has been chosen for simulations, 
with K = 150/c B T, l/7a = 2 nm and E min = -4k B T. 
With these parameter values, a good qualitative agree- 
ment is obtained, as displayed in figure [3] Cluster-size 
distributions are bimodal in both cases, with a monomer 
peak and clusters containing a few proteins (compare 
both our green curve in figure [7] and figure 4 of [70]). 
This good qualitative agreement supports our numerical 
approach. Similar cluster sizes {k « 10) have been ob- 



served at lower protein density but after photoactivation 
ofBR[7l]. 



4.2. Role of three-body elastic interactions 

It has been proved [3H H5] that the elastic energy for 
isotropic inclusions is dominated at lon g ra nge by the 
sum of the two-body terms of equation (3.2) considered 
so far and the three-body ones of equation (3.3 ). Here we 



take these three-body terms into account in simulations. 
Contrary to two-body forces, we did not include any cut- 
off at short distances, the hard-core repulsion playing this 
role. Note that in this paragraph, E m i n still refers to th e 
two-body energy at contact of U^" {2a) [equation (3.2)]. 



Simulating three-body forces consists of an important 
increase of the computational complexity, which grows 
from 0(N 2 ) to 0{N 3 ) operations per Monte Carlo step. 
For this reason, systems with N > 100 can hardly be 
simulated and we could perform a few simulations only, 
in order to check that cluster phases still exist in this case 
and thus that two-body repulsion dominates three-body 
forces, which can be attractive or repulsive (see also the 
Supporting Information) . 
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FIG. 4: Distributions of cluster sizes P(k) in the case of alge- 
braic long-range interaction including three-body interaction 
terms of equation ( |3.3[ ), for different parameter sets. In all 
cases, K = 30fcBT- The dot-dashed line has slope -2 for com- 
parison. Log-log coordinates. 



The principal features observed in the previous subsec- 
tion are conserved. Cluster phases appear to be stable 
at equilibrium and their distributions display a power- 
law-like regime until a maximum size. Examples are pro- 
vided in figure [4] When increasing the density, the mean 
cluster-size (fc) grows, even if we do not have enough 
points at our disposal to conclude that a proportional- 
ity regime still exists. Lowering E mul at fixed K also 
increases (fc). 

However, some differences are also noticeable. First, 
higher values of the energy at contact -Emm are required 
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to obtain distributions comparable to those of the pre- 
vious paragraph. Their are close to — 1.5/cbT instead of 
— Ak-gT. This indicates that three-body forces have an 
attractive contribution, on average, as already noticed 
in reference |43j and in the Supporting Information. In 
addition, at high <f), simulations show that a big cluster, 
containing a large fraction of the N particles, can tran- 
siently appear and subsequently disintegrate into smaller 
clusters. As a consequence, one distribution displays a 
secondary peak at large cluster sizes in the figure, for 
<f> = 0.2, E min = — 1.5fc B T and N — 150 (whereas it is 
absent if N = 100). This might be the signature of a 
transition between a cluster phase and a liquid-gas co- 
existence phase when increasing </> |69j . To clarify this 
point, simulating larger systems will be required. But 
again, simulating three-body forces is expensive and our 
statistics are poorer than in the previous sections. 



4.3. Modulation of long-range potentials 

We now explore the consequences of protein diversity 
in terms of long-range forces, for both kinds of pair po- 
tentials ( bi-ex ponential, equation (3.1), and algebraic, 
equation ( |3.2| ). Examples of modulated potentials are 
displayed in figure [5j 
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FIG. 5: Algebraic pair potentials [/J" 5 ™ (three upper curves) 
and bi-exponential ones [/| xp (lower curves) modulated at 
long range. This modulation takes into account the fact that 
cylindrical inclusions (in red) are up-down symmetric (an = 
0) while conical ones (in blue) are asymmetric (ub = 10°), 
as in the previous paragraphs. Here -E m i n = — 4&bT and 
l/7a = 2 nm for all curves, K — dOksT for the algebraic 
repulsion, E m&K — O.lksT, 0.05/cbT and for the three bi- 
exponential potentials, respectively (and l/j r = 16 nm). 



In the algebraic case, the modulation is ensured by 
taking a = (cylindrical inclusions thereafter) instead 
of a = 10° (conical inclusions) when calculating the con- 
stant C . This constant is essentially divided by 2 for large 
K values when a conical and a cylindrical inclusion inter- 
act. It becomes negative for two cylindrical inclusions, 



which feel a mutual attraction because of Casimir forces. 
As for bi-exponential potentials, we adapt both e a and e r 
so that .E m i n remains unchanged whereas E max is divided 
by 2 for cylindrical-conical interactions and E max = for 
cylindrical-cylindrical ones. In the latter case, we simply 
set e r = and adapt e a to set the value of Emm! The 
interaction is everywhere attractive. We denote by iV C yi. 
and iVcon. the numbers of cylindrical and conical parti- 
cles and ./V = iV cy i. + N con ,. We compare systems with 
a finite cylindrical- inclusion fraction x cy \. = N cy \jN to 
systems of N identical conical particles, as above. 

Figure [6] provides simulation snapshots. For small 
cylindrical-inclusion fractions x cy i, they are very simi- 
lar to the mono-particle case as studied above. Figure [7] 
shows examples of cluster-size distributions P(k) that are 
weakly affected by the presence of a small amount of 
cylindrical inclusions, at fixed <f>. Consequently, propor- 
tionality regimes are also preserved (figure [8j bottom). 

To interpret these numerical observations, we discuss 
them in the framework of the liquid-gas transition of bi- 
nary mixtures, even though in the present context of clus- 
ter phases, transitions are certainly not true thermody- 
namic transitions but crossovers [37|. We see a cluster 
as a liquid droplet co-existing with gas. In the liquid 
phase, conical and cylindrical inclusions are miscible be- 
cause they essentially feel the same short-range attrac- 
tion. This miscibility is visible on simulation snapshots 
(figure [6]). However, the attraction is slightly stronger 
for cylindrical inclusions. Their second virial coefficient 
B 2 = 7T f °° r[l — exp(— U(r))]dr is thus slightly weaker, 
and the transition occurs at a lower density. Conse- 
quently, clusters are bigger and enriched in cylindrical 
inclusions, as observed. The average fraction cc£ yl of 
cylindrical inclusions within clusters is > x cy \,. However, 
conical inclusions remain the majority in clusters and 
too large clusters are unstable, in spite of the cylindrical- 
inclusion fraction. The presence of a minority of such 
particles does not destabilize the cluster phase. 

However, we anticipate that there exists a value of x cy i. 
above which this mechanism will cease being valid. To 
test this hypothesis, we increased the fraction of cylin- 
drical inclusions up to x cy i — 0.4. Figure M shows our 
results for bi-exponential potentials, for both N = 100 
and 200 particles. As x cy \, grows, the cluster peak drifts 
to larger sizes in the bimodal distributions for N = 100. 
But the N = 200 simulations indicate that this large peak 
is in fact the superimposition of two peaks for x cy i. > 0.3. 
And the inspection of simulation snapshots indicates that 
there are indeed two kinds of clusters: mixed clusters as 
in the low a; cy i. regime, the size distribution of which 
matches the low x cy \, one; and one larger assembly that 
is highly enriched in cylindrical inclusions. This indi- 
cates the existence of a saturation phenomenon: beyond 
x cy i. — 0.3, additional cylindrical inclusions "precipitate" 
in a single large assembly with rare conical proteins. Note 
however that this large assembly is not stable: one ob- 
serves in simulations that this cluster can disintegrate 
and be reformed somewhere else. This instability seems 
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FIG. 6: Top: Snapshot of a cluster phase with 2V cy i. = 20 
cylindrical inclusions (in red) and N con . = 80 conical ones 
(in blue), in the case of bi-exponential pair potentials as in 
figure [5] The density is (f> = 1/50. Middle and bottom: 
Snapshots of cluster phases with AT cy i. = 10 cylindrical in- 
clusions and iV con . = 90 conical ones, in the case of the pair- 
wise algebraic repulsion with K = 150fcBT, -Emin = — 4£;bT 
and l/7a = 2 nm. Densities are 0=1/3 (middle) and 
4> — 1/7 (bottom). In all three cases, mono-particle systems 
(N C yi. = 0) give essentially similar snapshots. 
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FIG. 7: Distributions of cluster sizes P(k) in the case of the 
pairwise algebraic long-range repulsion with K — 150/cbT, 
E m i n — —AksT and l/j a = 2 nm, for two different densities 
<j> and two different cylindrical inclusion numbers N cy \. (N = 
100 particles in all cases). The green distribution is (weakly) 
bimodal. 



to be stronger at x cy \, = 0.35 than at x cy \. = 0.4. To 
sum up, if x cy i. is too large, many cylindrical inclusions 
aggregate in a single bulk liquid phase and remaining 
(cylindrical and conical) inclusions form a cluster phase 
as above, but with a lower effective x cy \.. 



5. MODULATION OF SHORT-RANGE 
POTENTIALS AND SPECIALIZATION OF 
CLUSTERS 

We now focus on the role of short-range potential mod- 
ulations in order to explore further the role of protein 
diversity. One of our major goals here is to validate 
on numerical grounds the analytical predictions of ref- 
erence [3 7) concerning protein cluster specialization: for 
a sufficiently large affinity difference between same- and 
distinct-family proteins, clusters become heterogeneous 
and contain essentially one protein family. We indeed 
consider systems with q families of proteins. Each family 
represents proteins, not necessarily identical, that have a 
preferential affinity E m i n at close range. There are typi- 
cally q = 10 3 families in a real membrane. The energy at 
contact of two particles in a same family i, called same- 
family particles hereafter, is denoted by S^^min < 0. 
The energy at contact of particles of two distinct fam- 
ilies i and j (distinct-family particles) is denoted by 

In this section, we focus on bi-exponential potentials 
[equation (3.1)] because large clusters are more abun- 



dant in this case due to their bimodal distribution. This 
provides better statistics when computing Binder cumu- 
lants below. We consider a single density, </> = 1/50, 
and the most symmetric case where the q families con- 
tain the same number M = N/q of particles. All same- 
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FIG. 8: Top: Distributions of cluster sizes in the case of bi- 
exponential pair interaction potentials as in figure [5] with 
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aicyi.- Comparison between N = 100 (soild lines) and TV = 
200 (dotted lines) for each value of x cy i.. Bottom: resulting 
mean cluster-size as a function of the density, exhibiting a 
proportionality regime. 



family energies at contact, Ei, 



are also identical, as 



well as distinct-family ones, E it j min . The energy bar- 
rier, i? max , is unchanged a nd e nergies at contact are 
set by tuning e a and e r in (3.1), as explained in Sec- 
tion |1 Following [37], we introduce a Flory parameter 



X - 6(£i 



min) = §AE |78J, illustrative of the 



difference of affinity between same- and distinct-family 
proteins. In this reference [37], it was demonstrated that 
there exist a threshold value \ c such that if \ < x C j 
particles of different families are mixed within clusters, 
whereas if x > X c \ they are demixed and clusters contain 
predominantly one family. Demixing is due to a com- 
petition between entropy of mixing and contact energy. 
Note that this transition (as well as the one at (f> c ) is 
certainly a simple crossover and not a true thermody- 
namical transition because of the finite cluster size |37j . 
The theoretical prediction is \ c — 2 for q = 2 and 

X c = 2^]n(q-l)^2lnq (5.5) 
q z 

for q > 3 (see also [75]). The great interest of equa- 



tion (5.5) lies in the fact that \ c grows like lnfc and AE 
remains of the order of /cbT even if q is large, making 
this scenario physically reasonable to account for cluster- 
specialization [37) . This is the reason why we wish to 
validate this result on numerical grounds. 
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FIG. 9: Snapshots of cluster phases with q = 2 families of 50 
particles each (blue and green; TV = 100), in the case of bi- 
exponential pair potentials with distinct-family energy at con- 
tact -Eij,mm = — 3.9&;bT (top, mixed clusters) and — 3.2&B7 1 



(bottom, demixed clusters), 
energy at contact is Ei t . 
<p= 1/50. 



In both cases, the same-family 
= — 4:ksT, and the density is 



Figure [9] shows two simulation snapshots in qualitative 
agreement with this prediction. We also simulated sys- 
tems with up to q = 5 families, showing the same agree- 
ment. To quantify further these qualitative observations, 
we determine the threshold value \ c by use of (modi- 
fied) Binder cumulants, an efficient tool to detect phase 
transitions numerically 73 ]. Our system being equiva- 
lent to a mean-field q- Potts model [37] , we naturally use 
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the following order parameter: 



m 



i 

3=1 



(5.6) 



limit). To ascertain that the measured thresholds were 
not corrupted by finite size effects, we also simulated a 
system of N = 210 particles for q = 3 (M = 70). We find 
j min = — 3.47/cb? 1 , in satisfactory agreement with the 
M = 50 result above. In addition, cluster size distribu- 



where xj is the fraction of family j proteins in a clus- tions for N — 150 or N = 210 are undistinguishable. 
ter and u) — exp(2i7r/g) is the g-th root of unity. The 
(modified) Binder cumulant is then defined as 



B, 



(m 4 ) fc 



(5.7) 



where averages are computed on clusters of size k. Fig- 
ure [10] displays the numerical values of B k for q = 2 as a 
function of k and Bij,min- All curves cross at Ef - min = 
— 3.56/cbT, which is the signature of a phase transi- 
tion for finite-size systems (the clusters in the present 
case) |73| . The corresponding value of the Flory parame- 
ter is x° = 2.64, which compares well to the theroretical 
value, x c — 2- 
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FIG. 10: Modified Binder cumulants Bk for q 
Same parameters as in figure [9] 



We followed the same procedure up to q = 5, by sim- 
ulating systems containing N — 50q particles (M = 50). 
The measured threshold values are listed in Table |TJ and 
compared to theoretical predictions. The agreement im- 



q 


Efj^min {k B T) 


x c 


X c , theory 


ratio 


2 


-3.56 


2.64 


2 


1.32 


3 


-3.44 


3.36 


2.77 


1.21 


4 


-3.36 


3.84 


3.30 


1.16 


5 


-3.27 


4.38 


3.70 


1.18 



TABLE I: Numerical threshold values of the energies at con- 
tact, -Ejj.min, and the Flory parameter \ as a function of the 
number of particle families q. Theoretical values of \ c come 
from equation (5.5 1. Here i?i,i, m in = — iksT, the system con- 



tains N = 50g particles, and the density is </> — 1/50. 



proves as q grows, thus supporting the theoretical ap- 
proach of reference |37| (in principle valid in the large q 
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FIG. 11: Cluster size distributions P(k). Same parameters as 
in figure [9] 



Cluster size distributions P(k) are qualitatively un- 
changed when tuning |-Ei Jimm |, as illustrated in figure 11 
(we get similar results for q > 3). However, when ch> 



creasing this parameter or the number of families q at 



fixed \Ei 



|, the typical cluster size decreases because 



the average short-range attraction decreases. In a sim- 
ilarly way, the fraction of monomers increases slowly 
when decreasing l-Eij.minl (with no apparent singularity 
at \Ef imin \). Note that cluster sizes used to compute 



Binder cumulants above (figure 10 ) were chosen near the 
maximum of P(k) to ensure significant statistical sam- 
pling. All these numerical observations confirm our pre- 
vious theoretical predictions [57| . 



6. DISCUSSION AND CONCLUSION 

Our first major goal was to show that cluster phases 
are robust with respect to the precise shape of interaction 
potentials, provided that they display both a short-range 
(~ 1 nm) attraction with a contact energy of a few kBT 
and a longer-range repulsion (~ 10 nm), with an energy 
barrier on the order of k^T or lower. We have focussed on 
two limiting potential shapes: a full, unscreened r -4 re- 
pulsion mediated by the elastic membrane; and a totally 
screened repulsion, decaying exponentially at long range. 
In the regime of parameters studied, clusters of size rang- 
ing from a few to several tens of particles co-exist with 
monomers above a critical protein density. We have also 
demonstrated that taking into account three-body forces 
propagated by the membrane, which have the same order 
of magnitude as two-body ones, does not destabilize clus- 
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ter phases. This suggests that cluster phases are generic 
in the cell membrane context. 

Our second core objective was to explore how protein 
diversity affects cluster phases. In this respect, we have 
modeled diversity by playing on interaction potentials, 
both at long and short range. At long range, we have 
proved that the introduction in the system of a minority 
(< 30 %) of up-down symmetric inclusions, experiencing 
mutual attraction instead of repulsion, does not destabi- 
lize cluster phases. At short range, we have confirmed 
previous analytical calcultations 37J: modulating inter- 
actions favors the segregation of proteins of different fam- 
ilies in distinct clusters, provided that proteins of a same 
family have a slightly higher affinity at contact than pro- 
teins of distinct families. Accordingly, the high protein 
concentration of a membrane (~ 50 % in mass) favors the 
co-localization of proteins of a same family in same clus- 
ters. Were the concentration be too low, proteins would 
essentially behave as a gas of monomers diffusing freely 
on the membrane. It would take a long time for them to 
encounter their partners in the membrane by diffusion. 
Clusterization and crowding due to high concentration 
thus favors interactions between them, which should be 
of biological relevance in terms, e.g., of faithful response 
to external stimuli [12] or of cell adhesion |14) . 

The present work intends to add some original contri- 
butions to a more ambitious program: deciphering cell 
membrane organization on physical grounds. However, 
even though we have clarified several points, we have 
been led to make some assumptions, some of which we 
intend to discuss now (see also the Supporting Informa- 
tion). The objective of the present work is not to be a 
definitive answer to the question of membrane organiza- 
tion, but rather to lay the foundations of an emerging 
scenario. 

We have explained why the effective forces propagated 
by the elastic membrane are dominated at large separa- 
tions by two- and three-body forces. But sub-dominant, 
higher-order forces exist that might play a role at in- 
termediate distances and thus influence cluster stability. 
However, we conjecture that cluster phases will remain 
stable in this case, because we have shown that the long- 
range elastic interaction remain repulsive at large sep- 
aration. Accordingly, distant clusters will feel a repul- 
sive force that will prevent their further coalescence in 
larger clusters. To definitely clarify this point, it would 
be necessary to fully take into account many-body con- 
tributions, which requires inverting a 3N x 3N matrix 
and to calculate its determinant at each Monte Carlo 
step [441 146) . For moderate system sizes, intensive sim- 
ulations on parallel computers might be able to address 



this question in a near future. 

Even though a live cell is out-of-equilibrium due to en- 
ergy consumption, experiments done in vitro after mem- 
brane removal from cells still reveal clusters [T2]. Since 
in this context, as in the one of artificial vesicles [15) . 
no active, energy-consumming processes can limit clus- 
ter size, the existence of finite clusters at equilibrium 
has to be understood. The mechanism explored in this 
study partakes of this approach. However, cluster phases 
properties, notably cluster-size distributions, are likely to 
be modified by active processes. Endocytosis/exocytosis 
and membrane recycling are likely to enhance cluster- 
ization because they fragment large assemblies of pro- 
teins. Furthermore, protein conformations change when 
they are activated [7U [75] , thus modifying their coupling 
with the membrane and consequently their interaction 
parameters. For example, photo-activation of BR in re- 
constituted proteoliposomes has been shown to promote 
the formation of clusters [7T]. Such ingredients could be 
incorporated in our simulations, e.g. by modeling pro- 
teins as systems with two internal states, the interac- 
tion parameters of which depend on their internal state. 
The perturbation of equilibrium properties could then be 
quantified. 

To clarify these points, our predictions need to be val- 
idated at the experimental level. As already pointed 
out [37], Forster Resonance Energy Transfer (FRET) [76] 
is adapted to quantify the typical distance between flu- 
orescently labeled proteins, and thus their degree of ag- 
gregation |15L I61j . By modulating physical parameters 
such as the inclusion asymmetry, the membrane rigidity 
or the hydrophobic mismatch [381 55 , it would in prin- 
ciple be possible to probe this degree of aggregation in 
model membranes such as vesicles fTS] or stacked sup- 
ported membranes [77]. Fluorescence Correlation Spec- 
troscopy (FCS) 76} should have similar capabilities, since 
it quantifies the density of entities diffusing in a mem- 
brane, a cluster essentially counting as a single entity. 



Supporting Information available 

I. Monte Carlo procedure; II. Poly-dispersity of pro- 
tein diameters and contact angles; III. On large cluster 
stability; IV. Comparing two- and three-body interac- 
tions for asymmetric inclusions; V. Rotational averag- 
ing of elastic interactions for anisotropic inclusions. This 
material is available free of charge via the Internet at 
http : //pubs . acs . org 
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